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Abstract 

We decribe a general result on the approximation of fixed points of 
operators between normed vector spaces allowing an explicit estimation 
of the error. This result is particularly suited for the approximation of 
invariant measures in dynamical systems and in particular by the Ulam 
method. 

We apply this result to implement an algorithm for the rigorous com- 
putation of invariant densities of piecewise expanding maps up to some 
error in the distance. 

We show how several related computational and numerical issues can 
be solved and show some computer experiment calculating rigorous ap- 
proximation of invariant measure and entropy of some one dimensional 
maps. 

1 Introduction 

Overview As it is well known, dynamical systems, are often very good models 
of natural phenomena. In many cases the presence of chaos is a problem for the 
simulation and the predictive understanding of the behavior of the system. On 
the other hand it is also known that the statistical behavior of the evolution of 
a system is quite stable and often from this point of view the evolution of the 
system obeys to some precise law which can be predicted and obtained from the 
description of the system. 

Several important features of the statistical behavior of a given system are 
"encoded" in the so called Physical Invariant Measure. Having quantitative 
information on the invariant measure can give information on the statistical 
behavior for the long time evolution of the system. Physical invariant measures 
are the ones which (in some sense that will be precised below) represent the 
statistical behavior of a large set of initial conditions. 

The problem of the existence and properties of such invariant measures has 
become a central area of research in the modern theory of Dynamical Systems. 
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A big part of the result is abstract and gives no quantitative precise information 
on the measure. This is of course a significant hmitation in the applications. 

What is said above strongly motivates the search for algorithms which are 
able to compute quantitative information on the physical measure. The prob- 
lem of approximating invariant measure of dynamical systems was quite widely 
studied in the literature, some algorithm are proved to converge to the real in- 
variant measure (up to errors in some given metrics) in some classes of systems, 
but results giving a rigorous bound on the error are relatively few. 

In this paper we will describe an algorithm which is able to approximate 
interesting invariant measures with a precise bound on the error of the approxi- 
mation and its practical implementation. The main theoretical ideas behind are 
quite general and in some sense simplify some previous approach (|18j). The 
practical implementation of the method and the necessary precise estimates 
are described here for the class of piecewise expanding maps. We also present 
some real computer experiment performing the rigorous computation on interval 
maps, and our solution to the nontrivial computational/numeric issues arising. 

Invariant measure and statistical properties Let X be a metric space, 
T : X X a. Borel measurable map and /i a T-invariant Borel probability 
measure. An invariant measure is a Borel probability measure fi on X such 
that for each measurable set A it holds n{A) = fi{T~^{A)). They represent 
equilibrium states, in the sense that probabilities of events do not change in 
time. 

A set A is called T-invariant if ^ A (mod 0). The system {X,T,fj.) 

is said to be ergodic if each T-invariant set has total or null measure. In such 
systems the famous Birkhoff ergodic theorem says that time averages computed 
along ^-typical orbits coincides with space average with respect to /x. More 
precisely, for any / e L^{X, fi) it holds 



for almost each x, where Sl = f + f o T + . . . + f o T"-\ 

This shows that in an ergodic system, the statistical behavior of observables, 
under typical realizations of the system is given by the spatial average of the 
observable with respect to fi. 

We say that a point x belongs to the basin of an invariant measure /i if (fT|) 
holds at X for each bounded continuous /. In case X is a manifold (possibly 
with boundary), a physical measure is an invariant measure whose basin has 
positive Lebesgue measure (for more details and a general survey see |25j). 

The transfer operator Let us consider the space SPM{X) of signed Borel 
probability measures on X. A function T between metric spaces naturally in- 
duces a function L : SPM{X) — > SPM{X) which is linear and is called transfer 
operator (associated to T). Let us define L: if /i e SPM{X) then L[^] is such 
that 




(1) 
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m{A) = ti{T-\A)). 



Measures which are invariant for T are fixed points of L, hence the computation 
of invariant measures (and many other dynamical quantities) very often is done 
by computing the fixed points (or other spectral information) of this operator 
(restricted to a suitable Banach subspace) . The most applied and studied strat- 
egy is to find a finite dimensional approximation for L (restricted to a suitable 
function space) reducing the problem to the computation of the corresponding 
relevant eigenvectors of a finite matrix. 

An example of this is done by discretizing the space X by a partition and 
replacing the system by a (finite state) Markov Chain representing the dynamics 
between sets of the partition (see precise definition in section [3]). Taking finer 
and finer partitions it is possible to obtain in some cases that the finite dimen- 
sional model will converge to the real one (and its natural invariant measure to 
the physical measure of the original system) . In some case there is an estimation 
for this speed of convergence (see eg. [5] for a discussion), but a rigorous bound 
on the error (and then a real rigorous computation) is known only in a few cases 
(piecewise expanding or expanding maps, see [3l 118)). 

Another approach is to consider a perturbation of the system by a small 
noise. The resulting transfer operator has a kernel and then can be approxi- 
mated by a finite dimensional one, by a kind of Faedo-Galerkin method and 
relevant eigenvectors are calculated (see e.g. [HIS]). 

Variations on the method of partitions are given in [BJlT], while in |23i a differ- 
ent method, fastly converging, based on periodic points is exploited for piecewise 
analytic Markov maps. In some interesting examples we can obtain the physical 
measure as limit of iterates of the Lebesgue measure /i = lim„_>.oo l"^) ■ To 
obtain rigorous bounds on the error the main point is to explicitly estimate the 
speed of convergence to the limit. This sometimes can be done using techniques 
related to decay of correlations (dHHl]). General, abstract results on the com- 
putability of invariant measures are given in [T3] (see also [H]). It is worth to 
remark that in these papers are shown also some negative result. Indeed, there 
are examples of computabl^ systems without computable invariant measures. 

Plan of the paper In section [5] we show a general result regarding the ap- 
proximation of fixed points for linear operators between Banach spaces. In this 
result fixed points are approximated extracting and exploiting as much infor- 
mation as possible from the approximated operator. This general result can be 
applied to the Ulam approximation method. In Section |3] we show how this can 
be done and we show an algorithm for the approximation of invariant measures. 
In particular we perform all the required estimations in the case of piecewise 
expanding maps (with bounded derivative). 

In Section |4] we show how to implement the algorithm in practice. In par- 
ticular we have to show a way to rapidly compute the steady state of a large 

^ Computable, here means that the dynamics can be approximated at any accuracy by an 
algorithm, see e.g. 1131 for precise definition. 
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Markov chain up to a prescribed error. We also discuss several other com- 
putational and programming issues, explaining how we have implemented the 
algorithm to perform real rigorous computations on some example of piecewise 
expanding maps. 

In Section|n]we show the result of some experiments. Here we consider piece- 
wise expanding maps, and the invariant measure is computed up to an error of 
less than 1% with respect to the Li distance. As an application we show a rigor- 
ous estimation of the entropy (by the Lyapunov exponent) of such maps. These 
estimations can be used as a benchmark for the validation of statistical methods 
to compute entropy from time series. We remark that for the experimental val- 
idation of these methods to understand how they converge fast to the real value 
of the entropy an exact estimate for the value is needed. Usually some simple 
system is considered, where the entropy value can be computed analitically. We 
give a method which can produce such estimation on nontrivial systems, where, 
due to the lack of a Markov structure, the convergence of statistical, symbolic 
methods is slower (see [IS] e.g.). 
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2 A general result on the approximation of fixed 
points 

Let us consider a restriction of the transfer operator to an invariant normed 
subspace (often a Banach space of regular measures) B CSPM{X) and let us 
still denote it by L:B B. Suppose that is possible to approximate L in a 
suitable way by another operator Lg of which we can calculate fixed points and 
other properties. 

Our extent is to exploit as much as possible the information contained in Lg 
to approximate fixed points of L. Let us hence suppose that f,fs(zB are fixed 
points, respectively of L and Lg. 



Theorem 1. Suppose that: 
a) \\Lgf-Lf\\B<e 
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h)3N such that \\L^{fs-f)\\B<k\\fs-f\\B 
c) Ls is continuous on B; 3C s.t.\fg e HL^yHe < C||g||e- 
Then 

\\f5-f\\B<2e Yl 

Remark 2. In the following we show how the above items a),h),c) are natural, 
in the context of approximating a fixed point of the transfer operator: a) means 
that in some sense Lg is an approximation of L. About b), the required N will 
be calculated from a description of Ls exploiting the fact that, under natural as- 
sumptions Ls asymptotically contracts the space of zero average signed measures 
in B. Remark that b) also means that there is no "projection" of f on other fixed 
points of Ls than fs . c) will be obtained from the way Ls is defined. 

Proof, (of Theorem [1]) 

11/5 -/He < 

< \\Lfifs~f)\\s + \\Lff^L^f\\B 

< l\\fs-f\\B + \\L^f-L''f\\B 



(applying item b)). Hence 



but 



hence 



||/,--/||B<2||Lf/-L^/||B 



N 

Lf-L^ = Y^Lr\Ls-L)L' 
fc=i 



N 

{L^^L«)f = Y.^^'\Ls-L)L^-'f 

k=l 
N 

= Y.LrHLs-L)f 

k=l 



by item c), hence 



AT 

||(L^-Lf)/||e < ^C^-'=||(L.--L)/|b 
fe=i 

< e ^ C 

ie[0,7V-l] 



by item a), and then 



\\fs-f\\5<2e V C\ 



ie[O.N-l] 



□ 
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3 Estimation with Li norm and Ulam method 



Let us suppose now that X is a manifold with boundary. Let us describe better 
the Ulam's Discretization method. In this method the space X is discretized 
by a partition If, (with k elements) and the system is approximated by a (finite 
state) Markov Chain with transition probabilities 

P,,^m{T-\lj)rMi)/m{Ii) (2) 

(where m is the normalized Lebesgue measure on X) and defining a correspond- 
ing finite-dimensional operator {Ls depend on the whole chosen partition but 
simplifying we will indicate it with a parameter S related to the size of the el- 
ements of the partition) we remark that in this way, to Ls it corresponds a 
matrix Pk = (Pij) . 

We remark that Ls can be seen in the following way: let Fs be the cr— algebra 
associated to the partition Is, then: 

Ls{f)^E{L{E{f\Fs))\Fs), (3) 

(see also [TH], notes 9 and 10 for some more explanations). Taking finer and 
finer partitions, in certain systems including for example piecewise expanding 
one-dimensional maps, the finite dimensional model converges to the real one 
and its natural invariant measure to the physical measure of the original system, 
see e.g. [8ll3[18]. 

For the sake of simplicity we will suppose that all sets Ij have the same 
measure: ra{Ij) = i. This will simplify some notation. The general case can 
be treated similarly. 

We now apply Theorem [T] to a more concrete case. Suppose that: 

• Ls is the Ulam approximation of L as defined above. 

• B^L^iX), 

• L satisfies a Lasota Yorke inequality of the type 

\\L^9\\b' <Xl\9\\B'+B\\g\\L,, (4) 

where B' is another Banach subspace of Li and A < 1. Inequalities of this 
type hold for several systems having some form of expansiveness or hyperbolicity 
(see [linilinillS and Theorem[H]e.g.). 

Remark 3. A useful remark is that the L- Y inequality gives us an upper bound 
for the B' norm of a fixed point. Indeed, if f is a fixed point, we have that 

\\f\\B'^\\LV\\B'<\l\f\\B'+B\\f\\r^,. 

Now, supposing — 1 and letting n go to infinity, we have that 

\\f\\B'<B. 
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Let us outline the main points which ahow the application of Theorem [T] in 
this case to construct an algorithm to approximate invariant measures. 

Remark 4. To estimate \ \Lsf — LfW^-^ as required in item a) of Theorem{^we 
remark that 

\\Lsf - Lf\\L, < \\Ls-L\\B'^LA\f\\B'; 

(where ||.||b'-!-Li is the operator norm, as an operator B' — )• Li) hence we can 
write the conclusion of the theorem in this way: 

Wfs - /IIl, <2NY^ n\Ls - LWb'^lMW- 



Ahother way to estimate \\Lsf — L/||lj which is slightly more complicated 
but can give better estimations is explained in Subsection \3.3.4\ 
In this setting, under suitable assumptions: 

11 the norm *s estimated by the coefficients of the L-Y inequality, indeed 

as seen in Remark\^' 

\\f\\B'<B 

(see also Section \3.3.1\ below): 

12 under suitable assumptions \\Ls ~ -L||e'-^Li is estimated a priori by the 

method of approximation (see Section \3.3.2\ below): 

13 the integer N relative to item b) in Theorem\^can be estimated by the matrix 

Pk relative to Ls. In practice what we need to do is to calculate the first 
N such that ||-Pfc^|y||i < ^- Here is the 1-operator norm of Pk, and 
V is the space of vectors with zero mean (see Section \3. 1\ below). 

14 Since B = Li{X) and we consider the Ulam approximation then C — 1 (see 

Section\K3J\ below). 

Now let us discuss more precisely Item 73, which is central in this approach 
and whose discussion is general. We discuss the other Items in the Subsection 
13.31 with precise estimations related to a particular family of cases: the piecewise 
expanding maps. 

3.1 About item 13 

To compute N we consider F = {/i e B\fi{X) — 0} and ||L^|v||li^Li. Since 
f ~ fs G V , if we prove 

we imply Item b) of theorem [T] In the Ulam approximation Ls is a finite rank 
operator, hence, once we fix a basis this is given by a matrix. The natural basis 
{/i, to consider is the set of characteristic functions of the sets in the 

partition Is. If Is = {Ii, 7fe} then fi = Ij.; after the choice of this basis, the 
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set of linear combinations of such characteristic functions can be identified with 
R*^. By a small abuse of notation we ill also indicate by V the set of zero average 
vectors in R'". 

To determine N we have to consider the matrix \ v associated to the action 
of Ls on the space of zero mean vectors with respect to this basis and compute 
its operator norm ||Pfc|y||i wher^ 

llPfclvlli = sup |P(i;)|i. 

\v\i = l 

By Equation [3] the behavior of Ls and its relation with Pk is described by 

f V ^ V ^ f ^ Lsif) 

where /: K.'^' — > Li is the trivial identification of a vector in M'^ with a piecewise 
constant function given by the choice of the basis. This implies that V/ G Li 

wlsWl.^l, < mwi- 

Indeed if / G Li , ||E(/|F5)||li < ||/||li and / is trivially an isometry. 

Remark that if / /dm = 0, then / E{f\Fs)dm = and converse, and hence 

WLsWWl, < ||Pfe|/-i(y)||i. 
Trivially the matrix corresponding to is P^^ . And then 

Summarizing, we can have an estimation of ||L^|y ||ij^_j.ij^by calculating a 
matrix AI approximating Pk\i-^(v) and ||M^||i. 

The algorithm will hence calculate HM-' lli for each integer j > 0, computing 
iteratively from A'P~^, until it finds some j for which HM-^Hi < ^ and 
output this j as the N required in item b) of Theorem [TJ 

3.2 The algorithm 

We now present informally the general algorithm which arises from the previous 
considerations for the approximation of invariant measures. More details on the 
implementation of each step are given in the following subsections. 

Algorithm 5. The algorithm hence works as follows: 

1. Input the map and the partition. 

2. Compute the matrix Pj. and the corresponding approximated fixed point fg 
of Lg up to some required approximation ei 

3. Compute AL, an estimation for \ \Lsf — L/||lj up to some error £2 
^|.|i will denote the Li norm on R". 
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4- Compute N such that item b) of Theorem{l\is verified as described in item 
13 above 

5. If all computations end successfully, output fs . 

All was said before allows us to state the following 

Proposition 6. I^^{ fs) is an approximation of one invariant measure of L up 
to an error e given by: 

e < ei + 2N{AL + £3) 

in the Li norm. 

3.3 Details in the piecewise expanding case 

We now enter in more details, showing how the previously explained algorithm 
works in a concrete but nontrivial family of cases, where all the required com- 
putations and estimations can be done. 
Let 

IImII = sup|^(0')l 

this is related to bounded variatiorH: if fi has density / then it is easy to show 
that ll^ll < 2war(/) + 2|/U 

In this case X — [0, 1], B' = {/i, ||/i|| < 00}. The dynamics we will consider 
is defined by a map satisfying the following requirements: 

Definition 7. A nonsingular function T : ([0, — >■ ([0, is said piece- 

wise expanding if 

• There is a finite set of points di ~ 0, (i2, — 1 such that rj^^;. c;.^^) is 

and for each i, J j^iy dx < 00. 

• infj.gjo,i] I -Da;?" I > 2 on the set where it is defined. 

We remark that usually it is supposed inf^-g^Q iDajTj > 1. We can always 
suppose that the derivative is bigger than 2 by considering some iterate of T (of 
course the physical measure of the iterate is the same). 

We suppose that the map is computable, in the sense that we can compute 
the probabilities Pij defined in [2] up to any given accuracy. This is the case for 
example, if the map is given by branches which are given by analytic functions 
with computable coefficients. 

It is now well known that piecewise expanding maps have a finite set of 
ergodic absolutely continuous invariant measure with bounded variation density. 

Such densities are also fixed points of the (Perron Frobenius) operatoiQ L : 
L^[0,1] i^O, 1] defined by 

^ Recall that the variation of a function g is defined as var(g) = 

SUP{a:;)gPinito subdivisions of [0, 1] J2i<n Isi^i) " 9{^i+l)\ < ^ 

■'Note that this operator corresponds to the above defined transfer operator, but it acts on 
densities instead of measures. 
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We now explain how to face all the points raised in the concrete implementation 
of Algorithm [SJ 

3.3.1 About Item II 

In this section we obtain an explicit estimation of the coefficients of the Lasota 
Yorke inequality for piecewise expanding maps. We follow the approach of jT^ . 
trying to optimize the size of the constants. 

Theorem 8. IfT is piecewise expanding ad above and is a measure on [0, 1] 

\\m < ^IImII + ^77^--r^M(i) + Mlj^D- 

Proof. Remark that 

Ze{(diA, + i)\ie{l,...,n-l)} 

since LpL gives zero weight to the points di {Lfi is absolutely continuous). 

For each such Z define <^cont to be linear and such that 0^ = on dZ, then 
define ip^ = cj)— (ji^^nv ■^i extend it to [0, 1] by setting it to zero outside 
Z . This is a continuous function. Moreover for each x £ Z 

\(Pcont\oo < 



min((ij - di+i) 
Thus 



z 



now remark that, on Z , -02 oT = [ ^^^ y + '"^^^^Tyi" ,then 

|W)I<|EM(^)' XT-HZ))\ + \T.^^( '^^\T^)^ Xt-hz))\ 

■ 2'^'°° -Ml) 



mm{di - di+i] 

<HiHfi^y)\ + 2|0UMl7|^l) + — iT^^Mi)- 

T' (T'y mm(di-di+i) 

^^"'^ is not C^, but it can be approximated as well as wanted by a 

function ip^ such that [V'e ~ J2z( ^t''^ )\°° and ^([V'e ~ J2z( ^t''^ )\) are as small 
as wanted. Hence 

Im(( ^, ) )l < \\m I ^, loo < IImII I0loo 

J inf T' 
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and 



rpit 



IW)I < M\ — I0loo + 2|0UMI7;^I)+ . ,r , A ^) 

infT' [T'Y imn{di-di+i) 



infT' 



-M1) + 2M|^|) 



□ 



Remark 9. We remark that from 

\\m<^M\+ \ 

infT' min(di - di+i) 

it is easy to estract 



M(l) + 2/i 



{T'Y 



inf T' 



niin(di - di+i) 



rpti 



l\2 



iMli 



Where ~ sup|/i((/))| coincides with the Li norm for a density of ii. 

|0|oo=l 

Remark 10. From now on, the following notation is going to be used throughout 
the paper 

2 T" 
5'-— 7:^-^ + 27^ . (5) 



(T'Y 



This constant plays a central role in our treatment and is the biggest obstruction 
in getting good estimates for the rigorous error. 

We remark that once an inequality of the form 

\\Lg\\B' <2X\\g\\s' +B'\\g\\i. 
is estabhshed (with 2A < 1) then, iterating, we have 

\\L^g\\B' < 2X\\Lg\\B'+B'\\Lg\\i < 2A(2A||Lg||f;' + S'HffHi) + 
and thus 



|i"5ll6' <2"A"||Lg|b,+ 



1 



1 - 2A 



B'\\g\ 



obtaining the inequahty in the form required at |4] and the coefficient 

1 



B = 



1 - 2A 



-B' 



which is important to estimate ||/|| in our algorithm. 
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3.3.2 About item 12 

As outlined before, on tlie interval [0, 1] we consider a partition made of intervals 
having length S. As remarked in Item 12 we need an estimation on the quality 
of approximation by Ulam discretization (a similar estimation was given in [181 
Lemma 4.1]). 

Lemma 11. For piecewise expanding maps, if Ls is given by the Ulam dis- 
cretization as explained before, for every f in BV[0, 1] we have that 

\\Lf-Lsf\\L,<S{2\ + l + B')\\f\\. 
where B' is defined in Remark ] 1 (A 
Proof. First of all, we have that trivially 

mr(E(/|^)) < varif ), 
moreover, for f E BV[0,1], holds 

mm)-f\\L,<S-varif). 
Thus for all / in BV[Q, 1] wc have 

\\{L - Ls)f\\L, < mLW\Ts)\Ts)) - L(E(/|J-,))|U, + mfl^s) - 
which implies, from the remark above 

- Ls)f\\L, < S ■ var{LE{f\Ts)) + WWl^s) - fh^- 
Using now Theorem [5] we have that 
WiL - Ls)f\\L, < S2\\\Eif\:Fs)\\ + SB'\\Eif\:Fs)\\L, + mm) - (6) 

but, from the properties of E(/|J') stated above we have that, since ||/|| > 
varif) and||/||>||/|U, 

WiL - Ls)f\\L, < S2X\\f\\ + SB'WfW, + S ■ Varif) < ^ + 1 + 

□ 

This gives us the last estimate which is needed by Item 3 of algorithm [5l 

Remark 12. Working with partitions subdivinding the interval [0, 1] in k equal 
subintervals, and recalling Remark\^we that, when f is an invariant measure, 
the inequality takes the form 
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3.3.3 About item 14 

It is easy to see tliat if Lg is given by the Ulam method 

\\Lsf\\L, < II/IIl,; 

indeed ||L/|U, < and \\E{f\Fs)\\L, < as seen in Section O and 

Ls comes from the composition of such functions. 



3.3.4 A little complication for a better estimation in our case 

We remark that from Equation [Bl we can also obtain in the same way, recalling 
that II/IIl, = 1 

\\{L-Ls)f\\L,<S{2X + l)\\f\\+SB'. 
By this (see remark [T2|) 



k 

This inequality can also be used in our algorithm to improve the precision of the 
estimation in the case of Li approximation by the Ulam approximation scheme, 
obtaining: 

ll/-/,ll,.<™H^£±i^. 



4 Implementing the algorithm 

In this section we explain the details in the implementation of our algorithm in 
the piecewise expanding case and some related numerical issue. The main points 
are the computation of a rigorous approximation of the related Markov chain 
and a fast method to approximate rigorously its steady state. We include some 
implementation and numerical supplementary remarks, which can be skipped 
at a first reading. 



4.1 Computing the Ulam approximation 

To compute the matrix of the Ulam approximation, we have developed an algo- 
rithm that computes the entries of the matrix with a rigorous error bound 
e. The computed matrix will be called 11. Our algorithm computes each entry 
and the error associated to each entry; the maximum of all these errors is e. 

If we partition the interval [0, 1] in k intervals /i, . . . , 7^ the «, j-th entry of 
Pk is given by equation ([2|). 

To compute the (i,j)-th entry, we partition the interval li in m smaller 
intervals Ji, . . . , Jm- This m has to be chosen to be bigger than the derivative 
of T\j. to ensure that the images T(J;) are shorter than the target interval Ij. 

For each / = 1, . . . , m we look at the images T{Ji) and we test whether this 
images are contained in Ij. Clearly there are three possible situations: 
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1. T(Ji) is contained in Ij; 

2. T{Ji) is contained in the coniplement of Ij] 

3. T{Ji) overlaps Ij but it is not contained. 

Clearly, [T]and[2]define unambiguosly whether the interval belongs to T^^IjO 
li or not. In case [T] we have that Ji belongs to T^^Ij, while in case [5] it belongs 
to the complement. Since we are computing m{T^^Ij H U) in case [T] we take 
into account the measure of J;, while in case [2] we simply forget about J;. 



Case 1 
Case 2 
Case 3 



li 






Ij 


-hh 
li 


Jl 




\ V- 

T(JI) Ij 


-h- 
li 


h- 


T(JI) 


H F-K 

Ij 


-h- 
li 

1 


h- 

1 


[ V- 


\ V- 

T(JI) Ij 

III 1 



The tricky part is how to deal with case |3l we would like to understand 
how much of the mass of the interval is sent into Ij by T. To perform the 
computation with a rigorous bound on the error we developed a simple iterative 
algorithm, with a stopping criterion, that computes the coefhcients of the Ulam 
approximation keeping track of the error done. First we compute a matrix 11' 
which approximates P^, then we "normalize" 11' in a way that the resulting 
matrix 11 represents a Markov chain. If k is big enough we can assume we are 
in one of the two following cases: the function is monotone on J; or has one 
discontinuity point in J;. 

We first study the monotone case. Since Ji is an "undecided" interval and 
the function is monotone, we know that there exists a point a which partitions Ji 
into two subintervals; one of those two subintervals is sent into Ij by T and the 
other is sent into the complement of Ij . Clearly, this a may not be computable 
exactly; to solve this issue what we do is to use a sort of bisection principle. We 
divide Ji in m smaller intervals; some of these intervals are sent into Ij by T, 
some others are sent into the complement of Ij and one of them is going to be 
"undecided" . Now, we can divide this interval into m smaller intervals; some of 
these intervals are sent into Ij by T and so on. We iterate this process and, if 
s is the number of steps, at the s-th step the "undecided" interval has length 

1 

k ■ 

The algorithm stops when this length is shorter than a prescribed length, leaving 
a small undecided interval; the length of this "undecided" interval adds up to a 
number eu. 
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Numerical Remark 1. To compute the images of the endpoints of the intervals, 
we are working with a library that computes exactly rational numbers (the 
GiNaC library [5]). Since all our examples are, in their continuity intervals, 
rational maps and the intervals have rational endpoints, the computation of the 
matrix is rigorous. 

Numerical Remark 2. The matrix has a special form which mirrors the graph 
of the function, i.e., if the interval does not contain a discontinuity the non 
zero entries of the matrix are distributed on each row near the image of the 
mid-point of the interval. As we already said, the number of nonzeros on a 
row is bounded from above by two times the maximum of the derivative in the 
interval. This a priori knowledge of the structure of the matrix permits us to 
fasten the computation, computing only the intervals "near" the image of the 
midpoint of the interval. To be more precise, we compute all the 2S intervals 
which are near the one containing the image of the midpoint, where S is an 
integer that bounds from above the derivative of the function in [0, 1]. 

In a similar way we study the case of the interval crossing a discontinuity 
point: if the discontinuity point belongs to J/ we divide J/ in smaller intervals. 
Only one of these intervals will contain the discontinuity point and we divide 
it into even smaller intervals; we continue this process until the interval that 
crosses the discontinuity has length smaller than the prescribed length. As 
before, the length of this "undecided" interval adds up to a number Sij such 
that 

Numerical Remark 3. If the interval J/ contains a discontinuity point, we can, 
again, infer the distribution of nonzeros on the row by thinking about how the 
maps behave: the only intervals for which the entries are non zero are those 
near to the images of the endpoints of the interval Ji. 

The maximum of all the is really important for all our estimates: we are 
going to denote it by e. 

We recall that we denote the matrix containing the computed coefficients 
by n', to distinguish it from P/j, the actual matrix of the Ulam discretization. 
In general, this matrix is a sparse matrix, whose number of nonzero in the j-th 
row depends on the derivative of the function in the interval Ij . To be sure that 
n' is a good approximation of a stochastic matrix we compute the sum of the 
elements for each row, subtract this number to 1 and spread the result uniformly 
on each of the nonzero elements of the row obtaining a new "markovized" matrix 
n. This assures us that the matrix 11 represents a Markov chain. If we look at 
the graph associated to 11', i.e., the graph where each nonzero element 11^^- is 
a vertex and there is an directed edge between i and j, it is easily seen that 11 
represents a Markov chain with exactly the same graph. 

If s is the maximum of the errors — Pij\^ we have that for each row, if we 
denote by nnz the number of nonzero elements of the row, that the sum of the 
elements of the row differs from 1 by at most nnz • e. So, if we spread the result 
uniformly on each of the nonzero elements of the row we have a new matrix 11 
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such that 

\llij-Pij \ <2-e, 
therefore, the matrix H is such that 

||Pfc-n||i <2-k-s. 

Numerical Remark 4. The markovization process is done numericahy working 
with integers; this is done taking advantage of the representation of floating 
points in double precision floating point arithmetics. Each one of the coefiicients 
of the matrix is between and 1; we multiply each cocflicicnt of the matrix by 
2^° and take its integer part. Wc then subtract all these integers to 2^° obtaining 
a number z; we divide the intcgc;r z by the number of nonzeros and we sum the 
result on each of the nonzeros. The remainder of the division of z by the mimbcr 
of nonzeros is then added to the first nonzero of the row. Then, we divide every 
coeflicient by 2^^; this permits us to ensure that the matrix has row sum 1. 
This procedure implies that, in practice, our markovization algorithm produces 
a matrix H such that 

||Pfc-n||i <3-k-s. 

The matrix 11 is the matrix we are going to work with and the "markoviza- 
tion" process ensures that the biggest eigenvalue of 11 is 1. 

4.2 Computing rigorously the steady state vector and the 
error 

In this section we only consider the case of a transitive Markov chain; the results 
here work fine also with non recurrent Markov chains but if the Markov chain 
is reducible there are some issues that can be solved but are not treated here. 
For us, mostly, irreducibility means that there is only one invariant density. 

Numerical Remark 5. Remark that if the Markov chain is reducible, we can 
still decompose the Markov chain using tools from graph theory to identify 
the different connected components of the graph of the Markov chain; this is 
not implemented in our software since the examples to which we are going to 
apply it are all topologically transitive. This implies transitivity in the Markov 
chain approximating them. Indeed, let li and Ij be the interior of two intervals 
of the partition, since the map is topologically transitive and the derivative is 
bounded away from zero, there exists an TV^ such that T'^'^{Ii) H Ij ^ 0, and 
this intersection is a union of intervals, with nonzero measure. Therefore, if 
we call N the maximum of all these Nij the matrix Pj^ has strictly positive 
entries and therefore the matrix Pk represents an irrcdiiciblc Markov chain. By 
the Perron- Frobenius theorem this implies that the steady state of the Markov 
chain is unique. 

We want to find the left eigenvector associated with the eigenvalue 1 of the 
matrix 11; we know that 11 has largest eigenvalue 1 since it has row sum 1 and 
therefore it admits as right eigenvector the vector (1, . . . , 1)-' . Moreover, since 11 
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represents a transitive Markov chain, the eigenvalue 1 has algebraic multiplicity 
1. 

To compute this eigenvector we use the power iteration method; given any 
initial condition bo, if we set 

h+i = &i • n, 

we have that bi converges to the eigenvector we are looking for; we want to 
find a stopping criterion that permits us to bound the numerical error of this 
operation from above. 

To do this we use build an enclosure for the eigenvector using an idea from 
the proof of the Perron-Frobenius theorem. The main idea of the proof of 
the Perron-Frobenius theorem [S] Theorem 1.1] is that a Markov matrix A 
(aperiodic, irreducible) contracts the simplex A of the vectors v having 1-norm 
1. 

This simplex is given by the convex combinations of the vectors ei, . . . , Cfc of 
the base; after I iterations, the diameter of the simplex is going to be bounded 

by 

P'e,-A'e,||i = ||A'(e,-e,)||i, 

for i,j= 1, . . . ,k and i ^ j. 

Therefore, to build an enclosure for the eigenvector, we iterate the matrix 
on the vectors {e^ — ej} for i,j — 1, . . . , fc and i ^ j; in practice, thanks to the 
triangle inequality 

\\A'{e. - e,)||i < \\A\e, - e,)||i + P'(ei - e,)\\,, 

we can iterate only the vectors of the form {ei — ej} for j — 2^ . . . ,k. 

Remark 13. Please remark that {ei — ej} for j = 2, . . . , fc is a base for the 
space of average vectors; therefore, while we build the enclosure we can also 
compute the integer N used in Theorem [7] as explained in \S.ll Later, in this 
section, we explain some subtle points about the computation of N. 

We implemented the enclosure algorithm in the following way; the program 
get as an input a threshold and we iterate the vectors {ei — ej}, with 

j = 2, . . . ,n and look at their 1-norm. For each j, we denote by Ij the integer 
such that 

||A'^ (ei - ej)\ \ < Enum- 

Let I the maximum over j of all the Ij : we have that 
diani(A'(A)) < max | (e^ - e^)] |i 

i,j = l,...,n 

< P'^(ei-ej)||i-h||A''(ei-e,)||i < 2£„„„. 

Therefore, for any initial condition bo, iterating it (renormalizing it) I times 
we get a vector contained in A' (A), whose numerical error is enclosed by 2enum- 
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Numerical Remark 6. To take into account the numerical error of the row- 
column multiplication we used a really rough estimate, knowing that 

||fioat(Ai;) - Av\\i < k ■ e,nach, 

where Cmach is the roundoff error, and by a easy induction 

llfloat(A'w) - A'-V\\i <l-k- ernach- 

Indeed, to take into account the numerical error, for each j, we take an Ij such 
that 

11^'^ (ei - ej)||l +lj ■ k ■ e,nach < £num- 

Therefore, we computed the eigenvector with a rigorous bound on the nu- 
merical error; what is left is to compute the rigorous error. 
From Section [3. 3. 21 we have an estimate for the quantity 

\\Lf-Lsf\\L^ 

depending on constants that, once the map is known, we can compute by hand. 

The main issue that remains to be solved is the computation of the number 
of iterations N needed for the Ulam approximation Lg to contract to 1/2 the 
space of average vectors as explained in Section [01 

In some way, we already assessed this question while we were computing the 
iterations of the simplex; the vectors ei — ej, with j = 1, . . . , are a base for 
the space of average vectors, so, while computing rigorously the eigenvector, 
we can compute also the number of iterations. But we have to be careful since 
we do not know the matrix Pk of the Ulam approximation Lg explicitly but we 
know only its approximation 11. To compute the number of iterations needed 
by Ls to contract the space of average vectors to 1/2, we need to use a chain 
of inequalities to give an upper bound to the norm of for j an integer. 

Indeed (see Section [XT]) 

mww < m n^' +n^')kiii < m - mvWi + m^wwi. 

By how we constructed the matrix n we have that ||Pfc — n||i < 2-fc-e where 
k is the size of the partition, e is the maximum of the errors in the computation 
of the coefficient of the matrix and 

ll^;?'-n-'|y||i = ||^/^-'(p.-n)n-i|y||i 

<Ell^rVlli-ll^fe-nklli-||n-Vlli 

2 = 1 

< 2 • J • fc • £, 

since \\Pk - n|y||i < 2 ■ k ■ e, |v||i < 1 and ||n''|y||i < 1 for every j,h. 
Therefore 

\\Pl\v\\<2-j-k-e+\\W\v\\,. 
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Numerical Remark 7. Please note that, if the approximation error is big and the 
size of the discretization is big, then, after a short number of iterates, 2 ■ j ■ k ■ s 
may be so big that the whole computation is destroyed. Therefore, we have to 
be careful, while choosing the error allowed on each coefficient of the matrix, to 
choose it small enough that no problem arises in this step. 

We recall now the sources of error in our computation, to make clear the 
last step of our program 

1. the discretization error, coming from the Ulam Discretization of the trans- 
fer operator; 

2. the approximation error: since we cannot compute exactly the matrix Pk, 
we have to approximate it by computing a matrix 11; 

3. the numerical error in the computation of the eigenvector. 

The rigorous error for the invariant measure is computed from all these errors, 
using the algorithm we will explain below. 

Remark 14. Please note that in the experiments we have chosen the diameter 
of the enclosure to be small enough so that the numerical error is small and the 
biggest source of uncertainity is the discretization error. 

Now, if we denote by / the invariant measure of the system, by Vk the 
eigenvector of the Ulam approximation Pk, by Ve the eigenvector of 11 and by v 
the eigenvector computed using 11, we have 

||/-«||i < \\f -Vk\\i + \\vk~v,\\i + \\v,~i\\^. 

In the first part of the section we showed how to bound the numerical error 

\\Ve - V\\l, 

item 13] in our list. 

We recall that, since we have satisfied items II, 12, 13 and 14 we get that the 
discretization error (item[T|) || / — is, by our theoretical estimates,: 

\\f -Vk\\i<2-N\\Lf -LsfWL^. 

The last thing we need to compute to get our rigorous estimate is a bound 

for 

\\vk - Weill, 

the approximation error, item [2] in our list. To do this we use again Theorem 
m computing the number of iterates A^J^I needed for 11 to contract to 1/2 the 
space of average vectors; then 

\\vk-Ve\\i<2N,\\Pk-nMvk\\i <m,-k-e. 

^please note that, if e is small, = N is expected. In the program we compute the two 
values indipendcntly, even if in general < N. 
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Finally, we have that 

\\f~d\\i<2-N\\Lf-Lsf\\Li+m,-k-s + \\v,~ij\\i, 
or explicitly 

11/ -will <2-N- , -B + m.-k-e+Wv.-dWi. 

k 

Remark 15. Please note that, if we use the estimate in section \3.3.4\ the last 
inequality becomes 

2A-B A- B -\- B' 

11/ - will < 2N + 4iVe • fc • e + ||we - will- 

k 

5 Rigorous computation of the Lyapunov expo- 
nent 

One of the uses of a rigorously computed density is the rigorous computation 
of the Lyapunov exponent of a map. The Lyapunov exponent at a point x, 
denoted by \(x), of a piecewise expanding map is defined by 



n 

\{x)= hrn -5]log((r)'(x)); 



n— S-+00 n ■ 

i=0 



by Birkhoff ergodic theorem, we have that, relative to an invariant measure /i, 
for /x-a.e. x we have that 

\{x) = (\og{\T'\)d^i^\. 



Our algorithm permits us to compute the density of an invariant measure 
with a rigorous error bound; by Young's inequality we have that 



1 /.I 



\og{\T'{x)\)f{x)dx- / \og{\T'{x)\)v{x)dx 







< max(log|T'(x)|)||/-^||i, 

a;e[0,l] 



where by a small abuse of notation v is identified with by a piecewise constant 
function representing it, see Section 13.11 Please note that since we only know 
11/ — w||i, which is a global quantity, there are no naive ideas to have a better 
estimation. 

Therefore, to compute the Lyapunov exponent, the only thing we have to 
do is to compute with a (relatively) small numerical error the integral 

1 

log(|r'(a;)|)wda;. 

Since f {x) > 1 for all x, we have that, in each of the continuity intervals of 
T, log(|T'(a;)|) is continuous and at least C^. 
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To compute the integral we used the rectangle rule, interpolating the function 
by a piecewise constant function on the homogeneous partition on which we 
computed the Ulam approximation; in the intervals in which the function has 
a discontinuity point we have a bigger error, of the order of the size of the 
interval times the variation of the function. Since the partitions we used in our 
examples are really small, the numerical error coming from the method and the 
error arising from the discontinuities are really small compared to the rigorous 
error. 

Numerical Remark 8. To ensure that the numerical error is much smaller than 
the rigorous errour (so that, rounding up the rigorous error, we can simply ig- 
nore it), we used the library iRRAM ([22 ) to evaluate the functions which we 
integrate. This library implements "Turing-rigorous" computation, i.e., com- 
putes rigorously the output of a function up to a prescribed error. In particular, 
we used this library to compute the function in such a way that all the digits 
of the double representation are rigorously computed. If we evaluate the func- 
tion on the midpoints of the partition, we have that integrating with respect 
to the invariant measure is, essentially, the inner product between the vector 
containing the values of the density and the vector containing the values of the 
function. Since all the values are computed up to the machine precision, we 
know that the relative error of such an operation is bounded from above by 
k ■ emach which, with respect to our data, gives an absolute error of the order of 
lO^^'^. This is negligible w.r.t. the rigorous error, which is of the order of 10"'^. 

6 Numerical experiments 

In this section we show the output of some complete experiments we made, 
where we used the programs described above. 

The code is now in an hybrid state: the routines that generate the matrix are 
written using the BOOST Ublas library and can run on almost any computer, 
at the same time the main computational issue is the enclosure method for 
the certified computation of the eigenvector, which needs more computational 
power. In this step the number of matrix-vector products is proportional to the 
size of the partition: in our examples the size of the partition is of the order 
of 10^. This forced us to implement and run our programs in a parallel HPC 
enviroment, using the library PETSc and running them on the CINECA Cluster 
SP6. 

The code for the programs, the matrices and the outputs of the cluster are 
found in the directory 

http: / /poisson.phc.unipi.it/^nisoli/invmeasure/ 
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6.1 The Lanford map 



For our first numerical experiment we chose one of the maps which were inves- 
tigated in [H]. The map T : [0, 1] [0, 1] given by 

T : X ^ 2x + l^x{^ ^ x) (mod 1). 

What seems to be a good approximation of the invariant measure of the map 
is plotted in figure 1 of the cited article. Since this map does not comply with 
the hypothesis of our article, i.e. there are some points where 1 < |-Da;T| < 2 
we study the map := T o T. Clearly, the invariant measures for T and 
coincide. We write the explicit form of T^; if f{x) = 2x + ■^x{l — x) and 
g{x) = 2x + ■^x{l — x) — 1 we have 

' fifix)) <x<l - ^VT7+| 

I figix)) |-^^/T7 <x<| -^7l7-| 
^ g{g{x)) f-y^yT7-|<x<l. 

In figure [1] you can see a plot of this map. 

First of all, please note that in every component where the map is continuous, 
the map is a polynomial. So, we can use exact arithmetics to compute the matrix 
n. Please note that the discontinuity points are irrational; this is taken care as 
we explained in Section 14.11 

The experiment is made on a homogeneous partition of [0, 1] in 1048576 
intervals. 

In figure [5] you can see the plot of the invariant measure we obtain through 
our method. 

We can compute directly from the explicit form of the mapping that 
A = 4/9 and that B' is bounded from above by 12.48. Below, some of the data 
(input and outputs) of our algorithm. 



£ <3-10-" 
£„nm < 0.0001 
N 17 
Ne 18 
B 112.24 
/ 25 

With all this data and the estimate in section 13.3.41 we can compute a rigorous 
error bound for the computed invariant measure: 

II/- w||i < 0.0081. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 1: The second iterate of the Lanford map 



Using all these estimates, we can compute the Lyapunov exponent for the 
original Lanford map, which is 0.658±0.008. For its second iterate, i.e., the map 
we studied with our experiment, the computed Lyapunov exponent is 1.315 ± 
0.015 ~ 2- (0.658±0.008) as we could expect from the properties of the Lyapunov 
exponent. 

6.2 A piece wise expanding map 



In this paragraph we study a sort of Lorenz map with bounded derivatives, 
T : [0, 1] ^ [0, 1] given by 



whose graph is plotted in figure [31 

Please note that also in this case, in every component where the map is 
continuous, the map is a polynomial. So, we can use exact arithmetics to 
compute the matrix 11. 

The experiment is made on a homogeneous partition of [0, 1] in 1048576 
intervals. 




(7) 
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1.4 
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0.4 - 



0.2 - 



I ' 1 ' ' ' ' ' ' ' 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 2: The invariant measure for the Lanford map 

The invariant measure we obtain through our method is plotted in figure ID 
For this map B' < 8.9 and A = 1/3. Befow, some of the data (input and 
outputs) of our algorithm. 



e < 10-12 

enura < 0.0001 

TV 13 

TVe 14 

B 26.7 

I 22 

With all this data and the estimate in section [3.3.41 we can compute a rigorous 
error bound for the computed invariant measure: 

II/- will < 0.0013. 

Using all these estimates, we can compute the Lyapunov exponent for this 
piecewise expanding map, which is 1.327 ± 0.002. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 3: Map (O 

6.3 A map without the Markov property 

The map T : [0, 1] ^ [0, 1] given by 

17 

T{x) = —X mod 1 (8) 
o 

whose graph is plotted in figure[5l This map does not enjoy the Markov property: 
since (17/5)*^ is never an integer the orbit of 1 is dense. A finite partition {7^} 
of [0, 1] is said to have the Markov property if f{Ii) O Ij ^ 9 imphcs Ij C f{Ii)- 
clearly, the partition must contain at least an interval of the form [ao, 1]. This 
implies that the partition must contain an interval of the form [ai, /(I)], which 
in turn implies that it contains an interval of the form [02, /^(l)] and so on, for 
every k, it must contain an interval of the form [afc,/'^(l)] and since the orbit 
of 1 is dense, this partition is made by an infinite number of intervals. 

Please note that also in this case, in every component where the map is 
continuous, the map is a polynomial. So, we can use exact arithmetics to 
compute the matrix 11. 

The experiment is made on a homogeneous partition of [0, 1] in 1048576 
intervals. 

The invariant measure we obtain through our method is plotted in figure |6l 
For this map B' < 17 and A = 5/17. Below, some of the data (input and 
outputs) of our algorithm. 
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1.2 



0.8 



0.6 



0.4 



0.2 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 4: The invariant measure for map ([7]) 



e < 1.75 -10 

ermm < 0.0001 

N 13 

Ne 14 

B 41.29 

I 20 



-10 



With ah this data and the estimate in section 13.3.41 we can compute a rigorous 
error bound for the computed invariant measure: 

\\f-i\\i < 0.013. 

Since the map is piecewise hnear, we know exactly the Lyapunov exponent 
of the map, which is ln(17) — hi(5). 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 5: Map ([8]) 



6.4 Another example 

We study the map T : [0, 1] — > [0, 1] given by 

r fx 

whose graph is plotted in figure [T] This map is really similar to map ([8]), but 
it is nonlinear in the two intervals [5/17, 10/17] and [10/17, 15/17], where it is 
defined by two branches of a polynomial of degree two. 

Please note that also in this case, in every component where the map is 
continuous, the map is a polynomial. So, we can use exact arithmetics to 
compute the matrix 11. 

The experiment is made on a homogeneous partition of [0, 1] in 1048576 
intervals. 

The invariant measure we obtain through our method is plotted in figure 
[51 Please note that, near 0.337 and 0.403 there are two small "staircase steps" 
which are visible only zooming the graph. 

For this map B' < 17.61 and A = 1/3. Below, some of the data (input and 
outputs) of our algorithm. 



< a; < ^ 
<x< ^ 

17 - 17 

^ <X<1. 



(9) 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 6: The invariant measure for map ([8]) 



< 2.19-10"" 

< 0.0001 
14 
15 

52.82 
21 

With ah this data and the estimate in section 13.3.41 we can compute a rigorous 
error bound for the computed invariant measure: 

II/- will < 0.004. 

Using all these estimates, we can compute the Lyapunov exponent for this 
piecewise expanding map, which is 1.219 ± 0.006. 



^num 

N 

B 
I 



6.5 Some data on the computations 

In this section we would hke to give some informations about the computations 
itselves; while the code that computes the matrix is subject to revision, we are 
going to reuse the code that computes the eigenvector rigorously. 
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Figure 7: Map ^ 

At the moment, the code that computes the matrix a purely serial code, 
which uses the rational numbers type (cl_RA) of the CLN library, exposed 
through the numeric type of the GiNaC library [2 ; the computations ran with- 
out any problem on our desktop and laptop computers, in a relatively short 
amount of time (of the order of 10000 seconds on a single processor). This rou- 
tine is highly parallelizable but, at the moment, the main computational issue 
arise from the rigorous computation of the eigenvalue, as the following data 
shows. 

To compute rigorously the eigenvector the total number of matrix-vector 
multiplications needed is 19.353.649 for the Lanford example, 18.336.646 for 
map dll), 20.264.432 for map jH]) and 20.567.470 for map The Lanford 

experiment was run on 96 SMT cpu (2 cpus per core) and the other experiments 
were run on 64 SMT cpus. The total running time for the Lanford code was 
6.263e -I- 04 seconds, the total running time for the experiment for map ([7|) was 
6.229e + 04 seconds, for experiment ([5]) was 6.1086e -I- 04 and for map ^ was 
6.1229e-H04 seconds. 

As you can see, the rigorous computation of the eigenvector (we underline 
the fact that it is rigorous, since the computation of the eigenvector itself is a 
matter of few seconds), is really time consuming; this is due to the fact that 
we are taking the powers of a sparse matrix. The main issue with such an 
operation is that the powers of a sparse matrix could be, in general, non sparse. 



29 



1.4 




0.6 - 
0.4 - 

0.2 - 

I ' 1 ' ' ' ' ' ' ' 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 8: The invariant measure for map ([9]) 

This means that, after some relatively fast iterations, the vector that we are 
multiplying is dense. 

There are a lot of tecniques that are being developed to compute matrix 
powers; it seems like it is one of the most important fields of applications of 
GPGPU computing; in the near future we are going to experiment if our code 
runs faster on them. 

On the other side, one of the big issues may be the speed of convergence of 
the Ulam method; in the future we are going to investigate other discretization 
schemes: if the dimension of the matrices is much smaller this would imply 
faster computations. 
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